# MONTE CARLO DATA GENERATED IN MATLAB
# stored in .csv format in this directory:
# C:\Users\franzese\Documents\Work\WPDOCS\Spatial Methods\Path Dep\PolMeth 2010\RData\
# with filenames of formats:
# beh_tT_mcM = behavior vector for observation-period T in trial M
# net_tT_mcM = network matrix for observation-period T in trial M

## version 16 July 2010, 10:30pm, Rob last editor

# siena estimation loop

## SET WORKING DIRECTORY TO WHERE MC DATA FROM MATLAB STORED
setwd('C:/Users/franzese/Documents/Work/WPDOCS/Spatial Methods/Path Dep/PolMeth 2010/RData')
library(RSiena)

NumActors <- 30
NumTrials <- 5
NumPeriods <- 4

results <- c(0,0)

#--------------------------------------------------------------------#
#------------------------------ LOOPY LOOP --------------------------#
for(i in 1:NumTrials)
{

# load mc data created in matlab
	net1 <- read.csv(file=paste("net_t1_mc",i,".csv",sep=""),sep=",",header=F)
	net1 <- as.matrix(net1)
	nets <- "net1"
	beh <- mat.or.vec(NumActors,NumPeriods)
	beh1 <- read.csv(file=paste("beh_t1_mc",i,".csv",sep=""),sep=",",header=F)
	beh1 <- as.matrix(beh1)
	beh[,1] <- beh1
if (NumPeriods>1.5)
	{
	net2 <- read.csv(file=paste("net_t2_mc",i,".csv",sep=""),sep=",",header=F)
	net2 <- as.matrix(net2)
	nets <- paste(nets,"net2",sep=",")
	beh2 <- read.csv(file=paste("beh_t2_mc",i,".csv",sep=""),sep=",",header=F)
	beh2 <- as.matrix(beh2)
	beh[,2] <- beh2
	}
if (NumPeriods>2.5)
	{
	net3 <- read.csv(file=paste("net_t3_mc",i,".csv",sep=""),sep=",",header=F)
	net3 <- as.matrix(net3)
	nets <- paste(nets,"net3",sep=",")
	beh3 <- read.csv(file=paste("beh_t3_mc",i,".csv",sep=""),sep=",",header=F)
	beh3 <- as.matrix(beh3)
	beh[,3] <- beh3
	}
if (NumPeriods>3.5)
	{
	net4 <- read.csv(file=paste("net_t4_mc",i,".csv",sep=""),sep=",",header=F)
	net4 <- as.matrix(net4)
	nets <- paste(nets,"net4",sep=",")
	beh4 <- read.csv(file=paste("beh_t4_mc",i,".csv",sep=""),sep=",",header=F)
	beh4 <- as.matrix(beh4)
	beh[,4] <- beh4
	}

### COPY NEXT IF{}, UPDATING TO 5.5 & *6, ETC, UP TO THE NUMBER OF TIME-PERIODS OBSERVED
### CLUNKY, I KNOW, BUT COULDN'T GET R TO CYCLE THROUGH LEFT-HAND SIDE VARIABLE-NAMES

if (NumPeriods>4.5)
	{
	net5 <- read.csv(file=paste("net_t5_mc",i,".csv",sep=""),sep=",",header=F)
	net5 <- as.matrix(net5)
	nets <- paste(nets,"net5",sep=",")
	beh5 <- read.csv(file=paste("beh_t5_mc",i,".csv",sep=""),sep=",",header=F)
	beh5 <- as.matrix(beh5)
	beh[,5] <- beh5
	}

# ESTIMATE WITH RSIENA

### CHANGE c() ENTRY OF NEXT LINE UP TO THE NUMBER OF OBSERVED NETWORKS; CLUNKY, I KNOW, BUT...
network <- sienaNet(array(c(net1,net2,net3,net4),dim=c(NumActors,NumActors,NumPeriods)))
behavior <- sienaNet(beh, type="behavior")
mydata <- sienaDataCreate(network,behavior)
myeff <- getEffects(mydata)
#fix(myeff)
myeff[,9]<- FALSE
myeff$include[myeff$functionName%in%"Amount of network change in period 1"] <- TRUE
myeff$include[myeff$functionName%in%"Amount of network change in period 2"] <- TRUE
myeff$include[myeff$functionName%in%"Amount of network change in period 3"] <- TRUE

### NEED ADD ANOTHER OF NEXT LINES FOR EACH PERIOD BEYOND 3RD, ANALOGOUSLY TO ABOVE
# myeff$include[myeff$functionName%in%"Amount of network change in period 4"] <- TRUE

myeff$include[myeff$functionName%in%"Similarity on behavior"] <- TRUE

myeff$include[myeff$functionName%in%"Amount of behavioral change in period 1 on behavior"] <- TRUE
myeff$include[myeff$functionName%in%"Amount of behavioral change in period 2 on behavior"] <- TRUE
myeff$include[myeff$functionName%in%"Amount of behavioral change in period 3 on behavior"] <- TRUE

### NEED ADD ANOTHER OF NEXT LINES FOR EACH PERIOD BEYOND 3RD, ANALOGOUSLY TO ABOVE
# myeff$include[myeff$functionName%in%"Amount of behavioral change in period 4 on behavior"] <- TRUE

myeff$include[myeff$type%in%"eval"&myeff$functionName%in%"beh. behavior average similarity"] <- TRUE

#fix(myeff)
mymodel <- sienaModelCreate(useStdInits = TRUE, projname='network_coevolution')
ans1 <- siena07(mymodel, data=mydata, effects=myeff)
thetavec.i <- ans1$theta
covtheta.i <- ans1$covtheta
setheta.i <- sqrt(diag(covtheta.i))
results.i <- cbind(thetavec.i,setheta.i)
results <- rbind(results,results.i)
}
#------------------------------ END OF LOOP --------------------------#
#--------------------------------------------------------------------#

#SAVE THE RESULTS; CHANGE OUTPUT NAME TO REFLECT SET-UP THIS MC
write.csv(results,file="RSiena_N30_T4_Rb1_Rn1_MC5.csv",append=FALSE,quote=TRUE,eol = "\n", na = "NA", row.names=F)